Handout 8 R: Multiple linear regression

The data

download.file("http://www.openintro.org/stat/data/evals.RData", destfile = "evals.RData")
load("evals.RData")
str(evals)
## 'data.frame':    463 obs. of  21 variables:
##  $ score        : num  4.7 4.1 3.9 4.8 4.6 4.3 2.8 4.1 3.4 4.5 ...
##  $ rank         : Factor w/ 3 levels "teaching","tenure track",..: 2 2 2 2 3 3 3 3 3 3 ...
##  $ ethnicity    : Factor w/ 2 levels "minority","not minority": 1 1 1 1 2 2 2 2 2 2 ...
##  $ gender       : Factor w/ 2 levels "female","male": 1 1 1 1 2 2 2 2 2 1 ...
##  $ language     : Factor w/ 2 levels "english","non-english": 1 1 1 1 1 1 1 1 1 1 ...
##  $ age          : int  36 36 36 36 59 59 59 51 51 40 ...
##  $ cls_perc_eval: num  55.8 68.8 60.8 62.6 85 ...
##  $ cls_did_eval : int  24 86 76 77 17 35 39 55 111 40 ...
##  $ cls_students : int  43 125 125 123 20 40 44 55 195 46 ...
##  $ cls_level    : Factor w/ 2 levels "lower","upper": 2 2 2 2 2 2 2 2 2 2 ...
##  $ cls_profs    : Factor w/ 2 levels "multiple","single": 2 2 2 2 1 1 1 2 2 2 ...
##  $ cls_credits  : Factor w/ 2 levels "multi credit",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ bty_f1lower  : int  5 5 5 5 4 4 4 5 5 2 ...
##  $ bty_f1upper  : int  7 7 7 7 4 4 4 2 2 5 ...
##  $ bty_f2upper  : int  6 6 6 6 2 2 2 5 5 4 ...
##  $ bty_m1lower  : int  2 2 2 2 2 2 2 2 2 3 ...
##  $ bty_m1upper  : int  4 4 4 4 3 3 3 3 3 3 ...
##  $ bty_m2upper  : int  6 6 6 6 3 3 3 3 3 2 ...
##  $ bty_avg      : num  5 5 5 5 3 ...
##  $ pic_outfit   : Factor w/ 2 levels "formal","not formal": 2 2 2 2 2 2 2 2 2 2 ...
##  $ pic_color    : Factor w/ 2 levels "black&white",..: 2 2 2 2 2 2 2 2 2 2 ...

Exercise 1

Is this an observational study or an experiment? The original research question posed in the paper is whether beauty leads directly to the differences in course evaluations. Given the study design, is it possible to answer this question as it is phrased? If not, rephrase the question.

This is observational study since it is not randomization nor the variables were controlled. I would phase the question as “is physical appearance (beatuy) of the college professor correlated to differences in course evaluation?”. Other factors were not collected in the study so it will be difficult to answer the question as it is originally phrased.

Exercise 2

Describe the distribution of score. Is the distribution skewed? What does that tell you about how students rate courses? Is this what you expected to see? Why, or why not?

hist(evals$score, prob=TRUE)
curve(dnorm(x, mean=mean(evals$score), sd=sd(evals$score)), add=TRUE)

The population for area is left right-skew and unimodal. The mean is 4.175, median is 4.3, and the range is 2.3 to 5. The students rate course higher than expected. If the professors were rated using a likert scale (1-5), I would expect the mean to be around 3.

The data does not follow a normal distribution.

  • Mean: 4.175

  • Median: 4.3

  • SD: 0.544

Exercise 3

Excluding score, select two other variables and describe their relationship using an appropriate visualization (scatterplot, side-by-side boxplots, or mosaic plot).

evals2 <- select_if(evals, is.numeric)
corrplot(cor(evals2), method = "number")

# average beauty decrease as age increase
fit <- lm(evals$bty_avg ~evals$age)
plot(evals$bty_avg ~evals$age)
abline(fit)

# average beauty increase as bty_m1upper increase
fit <- lm(evals$bty_avg ~evals$bty_m1upper)
plot(evals$bty_avg ~evals$bty_m1upper)
abline(fit)

# average beauty increase as cls_perc_eval increase
fit <- lm(evals$bty_avg ~evals$cls_perc_eval)
plot(evals$bty_avg ~evals$cls_perc_eval)
abline(fit)

# average beauty increase as bty_m1upper increase
fit <- lm(evals$bty_avg ~evals$bty_m1upper)
plot(evals$bty_avg ~evals$bty_m1upper)
abline(fit)

# average beauty increase as bty_m2upper increase
fit <- lm(evals$bty_avg ~evals$bty_m2upper)
plot(evals$bty_avg ~evals$bty_m2upper)
abline(fit)

# average beauty increase as bty_f1lower increase
fit <- lm(evals$bty_avg ~evals$bty_f1lower)
plot(evals$bty_avg ~evals$bty_f1lower)
abline(fit)

boxplot(evals$bty_avg ~ evals$rank)

boxplot(evals$bty_avg ~ evals$ethnicity)

boxplot(evals$bty_avg ~ evals$gender)

boxplot(evals$bty_avg ~ evals$language)

boxplot(evals$bty_avg ~ evals$cls_level)

boxplot(evals$bty_avg ~ evals$cls_profs)

boxplot(evals$bty_avg ~ evals$cls_credits)

boxplot(evals$bty_avg ~ evals$pic_outfit)

boxplot(evals$bty_avg ~ evals$pic_color)

Simple linear

fit <- lm(evals$score ~ evals$bty_avg)
plot(evals$score ~ evals$bty_avg)
abline(fit)

Exercise 4

Replot the scatterplot, but this time use the function jitter() on the y- or the x-coordinate. What was misleading about the initial scatterplot?

The initial plot has many points overlapping. The plots using jitter provided a better visual plot to determine the distribution. We can identify areas of higher concentration that was not easy to easy in the original plot.

fit = lm(evals$score ~ evals$bty_avg)
plot(jitter(evals$score) ~ evals$bty_avg)
abline(fit)

plot(evals$score ~ jitter(evals$bty_avg))
abline(fit)

plot(jitter(evals$score) ~ jitter(evals$bty_avg))
abline(fit)

Exercise 5

Let’s see if the apparent trend in the plot is something more than natural variation. Fit a linear model called m_bty to predict average professor score by average beauty rating and add the line to your plot using abline(m_bty). Write out the equation for the linear model and interpret the slope. Is average beauty score a statistically significant predictor? Does it appear to be a practically significant predictor?

m_bty <- lm(evals$score ~ evals$bty_avg)
plot(jitter(evals$score) ~ jitter(evals$bty_avg))
abline(m_bty)

summary(m_bty)
## 
## Call:
## lm(formula = evals$score ~ evals$bty_avg)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9246 -0.3690  0.1420  0.3977  0.9309 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    3.88034    0.07614   50.96  < 2e-16 ***
## evals$bty_avg  0.06664    0.01629    4.09 5.08e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5348 on 461 degrees of freedom
## Multiple R-squared:  0.03502,    Adjusted R-squared:  0.03293 
## F-statistic: 16.73 on 1 and 461 DF,  p-value: 5.083e-05

The equation for the regression line is:

\[score = 0.0666 * btyavg + 3.8803 \]

Slope interpretation:

When the average beauty (btyavg) of the professor goes up by 1, we would expect the rating evaluation (score) to go up by 0.0666 when all other variables are held constant.

When the average beauty (btyavg) is 0, the rating evaluation (score) is 3.8803.

Average beauty (bty_avg) is a statistically significant predictor. There is a lot of variance around the line and the slope is small. Average beauty explains 3.2% of the variance.

Exercise 6

Use residual plots to evaluate whether the conditions of least squares regression are reasonable. Provide plots and comments for each one.

plot(fit)

Conditions for linear regression:

  • Linearity

The regression plot shows many plots scattered and not closed to regression line. There is slightly linear positive linearity.

  • Normal residuals

The qq plot shows evidence of non-normality. Slightly left skewed distribution.

  • Constant variability

There seems to have constant variability.

  • Independence of Observations

There is no evidence of dependence.

Multiple linear regression

Exercise 7

P-values and parameter estimates should only be trusted if the conditions for the regression are reasonable. Verify that the conditions for this model are reasonable using diagnostic plots.

The residuals are nearly normal and data points have constant variability around regression line. The observations are independent of each other. The conditions for this model using least squares regression are reasonable.

Exercise 8

Is bty_avg still a significant predictor of score? Has the addition of gender to the model changed the parameter estimate for bty_avg?

The addition of gender to the model did change the parameter estimate for bty_avg.

m_bty_gen <- lm(evals$score ~ evals$bty_avg +  evals$gender)
summary(m_bty_gen)
## 
## Call:
## lm(formula = evals$score ~ evals$bty_avg + evals$gender)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8305 -0.3625  0.1055  0.4213  0.9314 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3.74734    0.08466  44.266  < 2e-16 ***
## evals$bty_avg     0.07416    0.01625   4.563 6.48e-06 ***
## evals$gendermale  0.17239    0.05022   3.433 0.000652 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5287 on 460 degrees of freedom
## Multiple R-squared:  0.05912,    Adjusted R-squared:  0.05503 
## F-statistic: 14.45 on 2 and 460 DF,  p-value: 8.177e-07
multiLines(m_bty_gen)

The equation for the regression line is:

\[score = 0.0742 * btyavg + gendermale * 0.1724 + 3.7473 \]

When the average beauty score of the professor goes up by 1, we would expect the rating evaluation to go up by 0.0742 when all other variables are held constant.

When the average beauty score is 0 and gender is female, the rating evaluation (score) is 3.8803 when all other variables are held constant.

When the average beauty score is 0 and gender is female, the rating evaluation (score) is 4.0527 when all other variables are held constant.

Exercise 9

What is the equation of the line corresponding to males? (Hint: For males, the parameter estimate is multiplied by 1.) For two professors who received the same beauty rating, which gender tends to have the higher course evaluation score?

For two professors, one male and one female, who received the same beauty rating, the male tends to have the higher course evaluation score.

score <- function(btyavg = NULL, gendermale = NULL)
{
  y <- 0.0742 * btyavg + 0.1724 * gendermale + 3.7473 
  y <- round(y, 3)
  return(y)
}

# male with beauty average of 2
pred <- score(2, 1)
pred
## [1] 4.068
# female with beauty average of 2
pred <- score(2, 0)
pred
## [1] 3.896
# male with beauty average of 5
pred <- score(5, 1)
pred
## [1] 4.291
# female with beauty average of 5
pred <- score(5, 0)
pred
## [1] 4.118
# male with beauty average of 8
pred <- score(8, 1)
pred
## [1] 4.513
# female with beauty average of 8
pred <- score(8, 0)
pred
## [1] 4.341

The equation for the regression line is:

When gender = Male, the equation will evaluate to:

\[score = 0.0742 * btyavg + 0.1724 * gendermale + 3.7473 \]

When gender = Female, the equation will evaluate to:

\[score = 0.0742 * btyavg + 3.7473 \]

Slope interpretation:

When the average beauty score of the professor goes up by 1, we would expect the rating evaluation to go up by 0.0742 when all other variables are held constant.

When the professor is male, we would expect the rating evaluation to go up by 3.7473 when all other variables are held constant.

Exercise 10

Create a new model called m_bty_rank with gender removed and rank added in. How does R appear to handle categorical variables that have more than two levels? Note that the rank variable has three levels: teaching, tenure track, tenured.

m_bty_rank <- lm(score ~ bty_avg + rank, data = evals)
summary(m_bty_rank)
## 
## Call:
## lm(formula = score ~ bty_avg + rank, data = evals)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8713 -0.3642  0.1489  0.4103  0.9525 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3.98155    0.09078  43.860  < 2e-16 ***
## bty_avg           0.06783    0.01655   4.098 4.92e-05 ***
## ranktenure track -0.16070    0.07395  -2.173   0.0303 *  
## ranktenured      -0.12623    0.06266  -2.014   0.0445 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5328 on 459 degrees of freedom
## Multiple R-squared:  0.04652,    Adjusted R-squared:  0.04029 
## F-statistic: 7.465 on 3 and 459 DF,  p-value: 6.88e-05
multiLines(m_bty_rank)

The equation for the regression line is:

When gender = Male, the equation will evaluate to:

\[score = 0.0678 * btyavg - 0.1607 * ranktenure track - 0.1262 * ranktenured + 3.9815 \]

Slope interpretation:

When the average beauty score of the professor goes up by 1, we would expect the rating evaluation to go up by 0.0678 when all other variables are held constant.

When the professor rank is tenure track, we would expect the rating evaluation to decrease by 0.1607 when all other variables are held constant.

When the professor rank is tenured, we would expect the rating evaluation to decrease by 0.1262 when all other variables are held constant.

When the professor rank is teaching, we would expect the rating evaluation to neither decrease or increase when all other variables are held constant.

The search for the best model

Exercise 11

Which variable would you expect to have the highest p-value in this model? Why? Hint: Think about which variable would you expect to not have any association with the professor score.

I would expect the following variables to have high p-value: color, outfit, language, ethnicity

m_full <- lm(score ~ rank + ethnicity + gender + language + age + cls_perc_eval 
             + cls_students + cls_level + cls_profs + cls_credits + bty_avg 
             + pic_outfit + pic_color, data = evals)

Exercise 12

Check your suspicions from the previous exercise. Include the model output in your response.

summary(m_full)
## 
## Call:
## lm(formula = score ~ rank + ethnicity + gender + language + age + 
##     cls_perc_eval + cls_students + cls_level + cls_profs + cls_credits + 
##     bty_avg + pic_outfit + pic_color, data = evals)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.77397 -0.32432  0.09067  0.35183  0.95036 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            4.0952141  0.2905277  14.096  < 2e-16 ***
## ranktenure track      -0.1475932  0.0820671  -1.798  0.07278 .  
## ranktenured           -0.0973378  0.0663296  -1.467  0.14295    
## ethnicitynot minority  0.1234929  0.0786273   1.571  0.11698    
## gendermale             0.2109481  0.0518230   4.071 5.54e-05 ***
## languagenon-english   -0.2298112  0.1113754  -2.063  0.03965 *  
## age                   -0.0090072  0.0031359  -2.872  0.00427 ** 
## cls_perc_eval          0.0053272  0.0015393   3.461  0.00059 ***
## cls_students           0.0004546  0.0003774   1.205  0.22896    
## cls_levelupper         0.0605140  0.0575617   1.051  0.29369    
## cls_profssingle       -0.0146619  0.0519885  -0.282  0.77806    
## cls_creditsone credit  0.5020432  0.1159388   4.330 1.84e-05 ***
## bty_avg                0.0400333  0.0175064   2.287  0.02267 *  
## pic_outfitnot formal  -0.1126817  0.0738800  -1.525  0.12792    
## pic_colorcolor        -0.2172630  0.0715021  -3.039  0.00252 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.498 on 448 degrees of freedom
## Multiple R-squared:  0.1871, Adjusted R-squared:  0.1617 
## F-statistic: 7.366 on 14 and 448 DF,  p-value: 6.552e-14
plot(m_full)

Exercise 13

Interpret the coefficient associated with the ethnicity variable.

The ethnicity not_minority variable increases the score by 0.123493 when all other variables are held constant.

Professor who are not from minority ethnicity tend to score 0.123493 more when other variables are held constant. Because the p-value > 0.5 the variable it’s not statistically significant.

Exercise 14

Drop the variable with the highest p-value and re-fit the model. Did the coefficients and significance of the other explanatory variables change? (One of the things that makes multiple regression interesting is that coefficient estimates depend on the other variables that are included in the model.) If not, what does this say about whether or not the dropped variable was collinear with the other explanatory variables?

The cls_profs has the highest p-value and was removed. The significance and coefficient of the some of the explanatory variables changed. This indicates that the dropped variable, cls_profs, might be collinear with variables whose significance has increased.

# remove cls_profs
m1 <- lm(score ~ rank + ethnicity + gender + language + age + cls_perc_eval 
             + cls_students + cls_level  + cls_credits + bty_avg 
             + pic_outfit + pic_color, data = evals)
summary(m1)
## 
## Call:
## lm(formula = score ~ rank + ethnicity + gender + language + age + 
##     cls_perc_eval + cls_students + cls_level + cls_credits + 
##     bty_avg + pic_outfit + pic_color, data = evals)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7836 -0.3257  0.0859  0.3513  0.9551 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            4.0872523  0.2888562  14.150  < 2e-16 ***
## ranktenure track      -0.1476746  0.0819824  -1.801 0.072327 .  
## ranktenured           -0.0973829  0.0662614  -1.470 0.142349    
## ethnicitynot minority  0.1274458  0.0772887   1.649 0.099856 .  
## gendermale             0.2101231  0.0516873   4.065 5.66e-05 ***
## languagenon-english   -0.2282894  0.1111305  -2.054 0.040530 *  
## age                   -0.0089992  0.0031326  -2.873 0.004262 ** 
## cls_perc_eval          0.0052888  0.0015317   3.453 0.000607 ***
## cls_students           0.0004687  0.0003737   1.254 0.210384    
## cls_levelupper         0.0606374  0.0575010   1.055 0.292200    
## cls_creditsone credit  0.5061196  0.1149163   4.404 1.33e-05 ***
## bty_avg                0.0398629  0.0174780   2.281 0.023032 *  
## pic_outfitnot formal  -0.1083227  0.0721711  -1.501 0.134080    
## pic_colorcolor        -0.2190527  0.0711469  -3.079 0.002205 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4974 on 449 degrees of freedom
## Multiple R-squared:  0.187,  Adjusted R-squared:  0.1634 
## F-statistic: 7.943 on 13 and 449 DF,  p-value: 2.336e-14
anova(m_full, m1)
## Analysis of Variance Table
## 
## Model 1: score ~ rank + ethnicity + gender + language + age + cls_perc_eval + 
##     cls_students + cls_level + cls_profs + cls_credits + bty_avg + 
##     pic_outfit + pic_color
## Model 2: score ~ rank + ethnicity + gender + language + age + cls_perc_eval + 
##     cls_students + cls_level + cls_credits + bty_avg + pic_outfit + 
##     pic_color
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    448 111.08                           
## 2    449 111.11 -1 -0.019722 0.0795 0.7781

Exercise 15

Using backward-selection and p-value as the selection criterion, determine the best model. You do not need to show all steps in your answer, just the output for the final model. Also, write out the linear model for predicting score based on the final model you settle on.

# remove cls_profs
m1 <- lm(score ~ rank + ethnicity + gender + language + age + cls_perc_eval 
             + cls_students + cls_level  + cls_credits + bty_avg 
             + pic_outfit + pic_color, data = evals)
summary(m1)
## 
## Call:
## lm(formula = score ~ rank + ethnicity + gender + language + age + 
##     cls_perc_eval + cls_students + cls_level + cls_credits + 
##     bty_avg + pic_outfit + pic_color, data = evals)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7836 -0.3257  0.0859  0.3513  0.9551 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            4.0872523  0.2888562  14.150  < 2e-16 ***
## ranktenure track      -0.1476746  0.0819824  -1.801 0.072327 .  
## ranktenured           -0.0973829  0.0662614  -1.470 0.142349    
## ethnicitynot minority  0.1274458  0.0772887   1.649 0.099856 .  
## gendermale             0.2101231  0.0516873   4.065 5.66e-05 ***
## languagenon-english   -0.2282894  0.1111305  -2.054 0.040530 *  
## age                   -0.0089992  0.0031326  -2.873 0.004262 ** 
## cls_perc_eval          0.0052888  0.0015317   3.453 0.000607 ***
## cls_students           0.0004687  0.0003737   1.254 0.210384    
## cls_levelupper         0.0606374  0.0575010   1.055 0.292200    
## cls_creditsone credit  0.5061196  0.1149163   4.404 1.33e-05 ***
## bty_avg                0.0398629  0.0174780   2.281 0.023032 *  
## pic_outfitnot formal  -0.1083227  0.0721711  -1.501 0.134080    
## pic_colorcolor        -0.2190527  0.0711469  -3.079 0.002205 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4974 on 449 degrees of freedom
## Multiple R-squared:  0.187,  Adjusted R-squared:  0.1634 
## F-statistic: 7.943 on 13 and 449 DF,  p-value: 2.336e-14
# remove cls_level
m2 <- lm(score ~ rank + ethnicity + gender + language + age + cls_perc_eval 
             + cls_students + cls_credits + bty_avg 
             + pic_outfit + pic_color, data = evals)
summary(m2)
## 
## Call:
## lm(formula = score ~ rank + ethnicity + gender + language + age + 
##     cls_perc_eval + cls_students + cls_credits + bty_avg + pic_outfit + 
##     pic_color, data = evals)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7761 -0.3187  0.0875  0.3547  0.9367 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            4.0856255  0.2888881  14.143  < 2e-16 ***
## ranktenure track      -0.1420696  0.0818201  -1.736 0.083184 .  
## ranktenured           -0.0895940  0.0658566  -1.360 0.174372    
## ethnicitynot minority  0.1424342  0.0759800   1.875 0.061491 .  
## gendermale             0.2037722  0.0513416   3.969 8.40e-05 ***
## languagenon-english   -0.2093185  0.1096785  -1.908 0.056966 .  
## age                   -0.0087287  0.0031224  -2.795 0.005404 ** 
## cls_perc_eval          0.0053545  0.0015306   3.498 0.000515 ***
## cls_students           0.0003573  0.0003585   0.997 0.319451    
## cls_creditsone credit  0.4733728  0.1106549   4.278 2.31e-05 ***
## bty_avg                0.0410340  0.0174449   2.352 0.019092 *  
## pic_outfitnot formal  -0.1172152  0.0716857  -1.635 0.102722    
## pic_colorcolor        -0.1973196  0.0681052  -2.897 0.003948 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4975 on 450 degrees of freedom
## Multiple R-squared:  0.185,  Adjusted R-squared:  0.1632 
## F-statistic:  8.51 on 12 and 450 DF,  p-value: 1.275e-14
# remove cls_students
m3 <- lm(score ~ rank + ethnicity + gender + language + age + cls_perc_eval 
             + cls_credits + bty_avg + pic_outfit + pic_color, data = evals)
summary(m3)
## 
## Call:
## lm(formula = score ~ rank + ethnicity + gender + language + age + 
##     cls_perc_eval + cls_credits + bty_avg + pic_outfit + pic_color, 
##     data = evals)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.78424 -0.31397  0.09261  0.35904  0.92154 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            4.152893   0.280892  14.785  < 2e-16 ***
## ranktenure track      -0.142239   0.081819  -1.738 0.082814 .  
## ranktenured           -0.083092   0.065532  -1.268 0.205469    
## ethnicitynot minority  0.143509   0.075972   1.889 0.059535 .  
## gendermale             0.208080   0.051159   4.067 5.61e-05 ***
## languagenon-english   -0.222515   0.108876  -2.044 0.041558 *  
## age                   -0.009074   0.003103  -2.924 0.003629 ** 
## cls_perc_eval          0.004841   0.001441   3.359 0.000849 ***
## cls_creditsone credit  0.472669   0.110652   4.272 2.37e-05 ***
## bty_avg                0.043578   0.017257   2.525 0.011903 *  
## pic_outfitnot formal  -0.136594   0.068998  -1.980 0.048347 *  
## pic_colorcolor        -0.189905   0.067697  -2.805 0.005246 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4975 on 451 degrees of freedom
## Multiple R-squared:  0.1832, Adjusted R-squared:  0.1632 
## F-statistic: 9.193 on 11 and 451 DF,  p-value: 6.364e-15
# remove rank
m4 <- lm(score ~ ethnicity + gender + language + age + cls_perc_eval 
             + cls_credits + bty_avg + pic_outfit + pic_color, data = evals)
summary(m4)
## 
## Call:
## lm(formula = score ~ ethnicity + gender + language + age + cls_perc_eval + 
##     cls_credits + bty_avg + pic_outfit + pic_color, data = evals)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8455 -0.3221  0.1013  0.3745  0.9051 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            3.907030   0.244889  15.954  < 2e-16 ***
## ethnicitynot minority  0.163818   0.075158   2.180 0.029798 *  
## gendermale             0.202597   0.050102   4.044 6.18e-05 ***
## languagenon-english   -0.246683   0.106146  -2.324 0.020567 *  
## age                   -0.006925   0.002658  -2.606 0.009475 ** 
## cls_perc_eval          0.004942   0.001442   3.427 0.000666 ***
## cls_creditsone credit  0.517205   0.104141   4.966 9.68e-07 ***
## bty_avg                0.046732   0.017091   2.734 0.006497 ** 
## pic_outfitnot formal  -0.113939   0.067168  -1.696 0.090510 .  
## pic_colorcolor        -0.180870   0.067456  -2.681 0.007601 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4982 on 453 degrees of freedom
## Multiple R-squared:  0.1774, Adjusted R-squared:  0.161 
## F-statistic: 10.85 on 9 and 453 DF,  p-value: 2.441e-15
backwards <- step(m_full, direction = "backward", trace=TRUE ) 
## Start:  AIC=-630.9
## score ~ rank + ethnicity + gender + language + age + cls_perc_eval + 
##     cls_students + cls_level + cls_profs + cls_credits + bty_avg + 
##     pic_outfit + pic_color
## 
##                 Df Sum of Sq    RSS     AIC
## - cls_profs      1    0.0197 111.11 -632.82
## - cls_level      1    0.2740 111.36 -631.76
## - cls_students   1    0.3599 111.44 -631.40
## - rank           2    0.8930 111.98 -631.19
## <none>                       111.08 -630.90
## - pic_outfit     1    0.5768 111.66 -630.50
## - ethnicity      1    0.6117 111.70 -630.36
## - language       1    1.0557 112.14 -628.52
## - bty_avg        1    1.2967 112.38 -627.53
## - age            1    2.0456 113.13 -624.45
## - pic_color      1    2.2893 113.37 -623.46
## - cls_perc_eval  1    2.9698 114.06 -620.69
## - gender         1    4.1085 115.19 -616.09
## - cls_credits    1    4.6495 115.73 -613.92
## 
## Step:  AIC=-632.82
## score ~ rank + ethnicity + gender + language + age + cls_perc_eval + 
##     cls_students + cls_level + cls_credits + bty_avg + pic_outfit + 
##     pic_color
## 
##                 Df Sum of Sq    RSS     AIC
## - cls_level      1    0.2752 111.38 -633.67
## - cls_students   1    0.3893 111.49 -633.20
## - rank           2    0.8939 112.00 -633.11
## <none>                       111.11 -632.82
## - pic_outfit     1    0.5574 111.66 -632.50
## - ethnicity      1    0.6728 111.78 -632.02
## - language       1    1.0442 112.15 -630.49
## - bty_avg        1    1.2872 112.39 -629.49
## - age            1    2.0422 113.15 -626.39
## - pic_color      1    2.3457 113.45 -625.15
## - cls_perc_eval  1    2.9502 114.06 -622.69
## - gender         1    4.0895 115.19 -618.08
## - cls_credits    1    4.7999 115.90 -615.24
## 
## Step:  AIC=-633.67
## score ~ rank + ethnicity + gender + language + age + cls_perc_eval + 
##     cls_students + cls_credits + bty_avg + pic_outfit + pic_color
## 
##                 Df Sum of Sq    RSS     AIC
## - cls_students   1    0.2459 111.63 -634.65
## - rank           2    0.8140 112.19 -634.30
## <none>                       111.38 -633.67
## - pic_outfit     1    0.6618 112.04 -632.93
## - ethnicity      1    0.8698 112.25 -632.07
## - language       1    0.9015 112.28 -631.94
## - bty_avg        1    1.3694 112.75 -630.02
## - age            1    1.9342 113.31 -627.70
## - pic_color      1    2.0777 113.46 -627.12
## - cls_perc_eval  1    3.0290 114.41 -623.25
## - gender         1    3.8989 115.28 -619.74
## - cls_credits    1    4.5296 115.91 -617.22
## 
## Step:  AIC=-634.65
## score ~ rank + ethnicity + gender + language + age + cls_perc_eval + 
##     cls_credits + bty_avg + pic_outfit + pic_color
## 
##                 Df Sum of Sq    RSS     AIC
## - rank           2    0.7892 112.42 -635.39
## <none>                       111.63 -634.65
## - ethnicity      1    0.8832 112.51 -633.00
## - pic_outfit     1    0.9700 112.60 -632.65
## - language       1    1.0338 112.66 -632.38
## - bty_avg        1    1.5783 113.20 -630.15
## - pic_color      1    1.9477 113.57 -628.64
## - age            1    2.1163 113.74 -627.96
## - cls_perc_eval  1    2.7922 114.42 -625.21
## - gender         1    4.0945 115.72 -619.97
## - cls_credits    1    4.5163 116.14 -618.29
## 
## Step:  AIC=-635.39
## score ~ ethnicity + gender + language + age + cls_perc_eval + 
##     cls_credits + bty_avg + pic_outfit + pic_color
## 
##                 Df Sum of Sq    RSS     AIC
## <none>                       112.42 -635.39
## - pic_outfit     1    0.7141 113.13 -634.46
## - ethnicity      1    1.1790 113.59 -632.56
## - language       1    1.3403 113.75 -631.90
## - age            1    1.6847 114.10 -630.50
## - pic_color      1    1.7841 114.20 -630.10
## - bty_avg        1    1.8553 114.27 -629.81
## - cls_perc_eval  1    2.9147 115.33 -625.54
## - gender         1    4.0577 116.47 -620.97
## - cls_credits    1    6.1208 118.54 -612.84

The best model is

\[score = 3.90703 + 0.16382 * ethnicitynot minority + 0.20260 * gendermale - 0.24668 * languagenon-english - 0.00692 * age + 0.00494 * cls_perc_eval + 0.51721 * cls_creditsone credit + 0.04673 * bty_avg - 0.11394 pic_outfitnot formal - 0.18087 * pic_colorcolor\]

Exercise 16

Verify that the conditions for this model are reasonable using diagnostic plots.

The conditions for this model are reasonable given that the diagnostic plots shows that residuals are nearly normal.

plot(m4)

qqnorm(backwards$residuals)
qqline(backwards$residuals)

plot(backwards$residuals)
abline(h = 0, lty = 3)

Assuming variable independence, the conditions for normality and homoscedasticity are met since in the qq-plot the data points fall along the normal line and the spread of the residuals along the zero line does not show a pattern.

Exercise 17

The original paper describes how these data were gathered by taking a sample of professors from the University of Texas at Austin and including all courses that they have taught. Considering that each row represents a course, could this new information have an impact on any of the conditions of linear regression?

This new information would potentially affect the linear regression because there could be overlapping data ormulticollinearity. This might potentially violate the assumption of independence.

Exercise 18

Based on your final model, describe the characteristics of a professor and course at University of Texas at Austin that would be associated with a high evaluation score.

The characteristics of a professor and course at University of Texas at Austin that would be associated with a high evaluation score:

A male professor typically receive higher evaluation. A professor that is not part of ethic minority group typically receive higher evaluation. A professor with higher average beauty score typically receive higher evaluation. A professor that teaches a one-credit course typically receive higher evaluation. A professor that received an education in a school that taught in english typically receive higher evaluation. A younger professor typically receive higher evaluation. A professor that uses a photo that has color typically receive lower evaluation. A professor that does not wear formal outfi typically receive lower evaluation.

Exercise 19

Would you be comfortable generalizing your conclusions to apply to professors generally (at any university)? Why or why not?

I won’t be comfortable in generalizing the conclusion as the sample may not be representative of entire population